#updated 1/8/2015 to work with R 3.1.2 and lme4 1.1-7
# also commented out graphics and summaries
#updated 11 27.2010
# reanalysis using BIC model selection
#=======================================================
#READING FILES
#=======================================================
grow=read.csv("July 2007 census and biomass2 edited.csv",header = T,na.strings=c("NA","na"))
#str(grow)
#converting variables in "grow" data to numeric/categorical variables
#is.factor(grow$comp)
#table(grow$comp)
grow$comp = as.factor(grow$comp)
grow$herb = as.factor(grow$herb)
grow$species = as.factor(grow$species)
#grow$no.leaves07=as.numeric(grow$no.leaves07)
grow$veg.bio = as.numeric(grow$veg.bio)
grow$bolt08= as.factor(grow$bolt)

# figure out which rows have the ugly "" factor levels
#grow[grow$Comp=="",] # its alot, but they are usually in blank rows
grow[grow$Comp=="",5:10] = NA
# relevel competition factor
grow$Comp = factor(grow$Comp,levels=c("Low","Medium","High"))
grow$Species = factor(grow$Species,levels=c("BT","TT"))
grow$Herb = factor(grow$Herb,levels=c("Ambient","Reduced"))
# I think this next step is overkill - there are many NA's that are in columns we don't use
#grow3=na.omit(grow)  #removed rows with NAs. in this case, dead plants are removed from the 
#analysis t
# discard rows with missing treatment values - mostly empty rows
grow = grow[complete.cases(grow[,5:10]),]

# 
grow = grow[complete.cases(grow[,c(11,18,19,26)]),]
grow2=grow[grow$no.leaves08>0 & grow$no.leaves07>0,]


#======================================================
#Visualizing the data
#======================================================

##Separating the data by species (bull vs. tall)
grow.bull=grow2[grow2$species==0,]
grow.tall=grow2[grow2$species==1,]
grow.tall$log.no.leaves08=log(grow.tall$no.leaves08)
grow.tall$log.no.leaves07=log(grow.tall$no.leaves07)
grow.bull$log.no.leaves08=log(grow.bull$no.leaves08)
grow.bull$log.no.leaves07=log(grow.bull$no.leaves07)


# par(mfrow=c(2,2))
# hist(grow.bull$no.leaves07)
# hist(grow.tall$no.leaves07)
# hist(grow.bull$no.leaves08)
# hist(grow.tall$no.leaves08)
# #hist(grow.tall$no.leaves08[grow.tall$bolt08==1], breaks = 2000, xlim = c(0,100))
# hist(grow.tall$no.leaves08[grow.tall$bolt08==0],breaks=20)
# 
# #plotting the scatter plot
# windows()
# par(mfrow=c(1,2))
# plot(log(grow.bull$no.leaves08),log(grow.bull$veg.bio), main="bull thistle")
# plot(log(grow.tall$no.leaves08),log(grow.tall$veg.bio), main="tall thistle")
# plot(grow.bull$no.leaves08,grow.bull$veg.bio, main="bull thistle")
# plot(grow.tall$no.leaves08,grow.tall$veg.bio, main="tall thistle")

# ===========================================================================
# modelling biomass as a function of leaves - Bull Thistle
# ===========================================================================
summary(gb.lm <- lm(veg.bio~no.leaves08*Herb*Comp, data = grow.bull))
#plot(gb.lm) # serious heteroscedasticity
summary(lgb.lm <- lm(log(veg.bio)~no.leaves08*Herb*Comp, data = grow.bull))
# point 22 is a serious outlier - eliminate it
remove = which(rownames(grow.bull)==22)
summary(lgb.lm <- lm(log(veg.bio)~no.leaves08*Herb*Comp, data = grow.bull[-remove,]))
# still some signs of non-linearity - log the number of leaves as well
summary(lgb.lm2 <- lm(log(veg.bio)~log(no.leaves08)*Herb*Comp, data = grow.bull[-remove,]))
# OK that has removed all the problematic outliers and cleaned up the heteroscedasticity

# now do it using lme4 to include random effects and bolting
# USE BIC model selection
library(lme4) #works with lme4 v1.1-7 and R 3.1.2
# make a unique "plot" identifier
grow.bull$plot = paste(grow.bull$block,grow.bull$plot,sep=".")

lgb.mm0 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (1|block) + (1|plot),data=grow.bull[-remove,])
lgb.mm1 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (1|block),data=grow.bull[-remove,])
lgb.mm2 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (Herb|block) + (1|plot),data=grow.bull[-remove,])
lgb.mm3 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (1|block) + (Comp|plot),data=grow.bull[-remove,])
lgb.mm4 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (1|plot),data=grow.bull[-remove,])
sapply(list(lgb.mm0,lgb.mm1,lgb.mm2,lgb.mm3,lgb.mm4),BIC)
# BIC values to pick random effect model
# go with block effect only, check if it makes a difference to model selection later. 

# create model set
bb.models = list("~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + 	Herb:Comp + Herb:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + 	Herb:Comp + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + Herb:bolt08 + Comp:bolt08 ",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp + log(no.leaves08):bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Herb + log(no.leaves08):bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + Herb:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + Herb:Comp",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp + log(no.leaves08):bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Herb",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + Herb:Comp",
		"~log(no.leaves08) + Herb + Comp + bolt08 + Herb:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08",
		"~log(no.leaves08) + Herb + bolt08 +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + Herb:bolt08",
		"~log(no.leaves08) + Herb + Comp + log(no.leaves08):Comp +  log(no.leaves08):Herb + Herb:Comp",
		"~log(no.leaves08) + Comp + bolt08 + log(no.leaves08):Comp + log(no.leaves08):bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + bolt08 +  log(no.leaves08):Herb + log(no.leaves08):bolt08",
		"~log(no.leaves08) + Herb + Comp + log(no.leaves08):Comp +  log(no.leaves08):Herb",
		"~log(no.leaves08) + Comp + bolt08",
		"~log(no.leaves08) + Herb + bolt08",
		"~log(no.leaves08) + Herb + Comp",
		"~log(no.leaves08) + Comp + bolt08 + log(no.leaves08):Comp + log(no.leaves08):bolt08",
		"~log(no.leaves08) + Herb +  log(no.leaves08):Herb",
		"~log(no.leaves08) + Comp + log(no.leaves08):Comp",
		"~log(no.leaves08) + bolt08 + log(no.leaves08):bolt08",
		"~log(no.leaves08) + Herb",
		"~log(no.leaves08) + Comp",
		"~log(no.leaves08) + bolt08",
		"~log(no.leaves08)",
		"~1")



bb.fits = lapply(bb.models,function(ff)lmer(paste("log(veg.bio)",ff,"+ (1|block)"),data=grow.bull[-remove,],REML=FALSE))
bic = sapply(bb.fits,BIC) 
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(bb.models),dbic,wbic)[order(dbic),]

# make some descriptive plots 
# library(lattice)
# xyplot(log(veg.bio)~log(no.leaves08)|bolt08+Comp,data=grow.bull[-remove,],groups=Herb)
# xyplot(log(no.leaves08)~log(no.leaves07)|bolt08+Comp,data=grow.bull[-remove,],groups=Herb,
# 	panel=function(...){
# 		panel.xyplot(...)
# 		panel.abline(a=0,b=1)
# 	}
# )


# ===================================================================
# growth model using veg.bio.07 - Bull thistle
# ==================================================================
nd = grow.bull[-remove,]
nd$no.leaves08  = nd$no.leaves07
nd$bolt08 = factor(0,levels=levels(grow.bull$bolt08))
results = matrix(NA,nrow=nrow(grow.bull[-remove,]),ncol=length(bb.fits))
for (i in 1:length(bb.fits)){
	# generate predictions from a lmer model using the data from inside the model, mostly. 
	fm = bb.fits[[i]]
  results[,i] = predict(fm,newdata=nd)
}
# model averaged prediction for veg.bio07 parameter
grow.bull[-remove,"veg.bio07"] = apply(results,1,function(x)sum(x*wbic))
# xyplot(log(veg.bio)~veg.bio07|Herb+Comp,data=grow.bull[-remove,],groups=bolt08,
# 	panel=function(...){
# 		panel.xyplot(...)
# 		panel.abline(a=0,b=1)
# 	}
# )


# plot(log(veg.bio)~veg.bio07,data=grow.bull, xlab="Log (Biomass) in t = 1", ylab="Log (Biomass) in t = 2")
# points(log(veg.bio)~veg.bio07,data=grow.bull[grow.bull$bolt08==1,], col="red")
# 
# xyplot(log(veg.bio)~veg.bio07|Herb+Comp,data=grow.bull,groups=bolt08)
# xyplot(no.leaves08~no.leaves07|Herb+Comp,data=grow.bull[-remove,],groups=bolt08)
# 
# create model set
bg.models = list("~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + veg.bio07:bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + veg.bio07:bolt08 + 	Herb:Comp + Herb:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + veg.bio07:bolt08 + 	Herb:Comp + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + veg.bio07:bolt08 + Herb:bolt08 + Comp:bolt08 ",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp + veg.bio07:bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Herb + veg.bio07:bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 +  veg.bio07:Herb + veg.bio07:bolt08 + Herb:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp + veg.bio07:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + Herb:Comp",
		"~veg.bio07 + Herb + Comp + bolt08 + Herb:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08",
		"~veg.bio07 + Herb + bolt08 +  veg.bio07:Herb + veg.bio07:bolt08 + Herb:bolt08",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp +  veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Comp + bolt08 + veg.bio07:Comp + veg.bio07:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + bolt08 +  veg.bio07:Herb + veg.bio07:bolt08",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp +  veg.bio07:Herb",
		"~veg.bio07 + Comp + bolt08 + veg.bio07:Comp + veg.bio07:bolt08",
		"~veg.bio07 + Herb + bolt08",
		"~veg.bio07 + Herb + Comp",
		"~veg.bio07 + Comp + bolt08",
		"~veg.bio07 + Herb +  veg.bio07:Herb",
		"~veg.bio07 + Comp + veg.bio07:Comp",
		"~veg.bio07 + bolt08 + veg.bio07:bolt08",
		"~veg.bio07 + Herb",
		"~veg.bio07 + Comp",
		"~veg.bio07 + bolt08",
		"~veg.bio07",
		"~1")

grb.mm0 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(1|block) + (1|plot),data=grow.bull[-remove,],REML=TRUE)

grb.mm1 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(1|plot),data=grow.bull[-remove,],REML=TRUE)
grb.mm2 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(1|block),data=grow.bull[-remove,],REML=TRUE)
grb.mm3 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(Herb|block) + (1|plot),data=grow.bull[-remove,],REML=TRUE)
# this model seems to have convergence problems, but randomly deleting a point makes it go away, and the conclusions don't change
grb.mm4 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(1|block) + (Comp|plot),data=grow.bull[-remove,],REML=TRUE)

sapply(list(grb.mm0,grb.mm1,grb.mm2,grb.mm3,grb.mm4),BIC)
# OK, don't need the block random effect, grb.mm1 is it
bg.fits = lapply(bg.models,function(ff)lmer(paste("log(veg.bio)",ff,"+ (1|plot)"),data=grow.bull[-remove,],REML=FALSE))
bic = sapply(bg.fits,BIC)
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
loglik=sapply(bg.fits,logLik) #unlist(sapply(bg.fits,function(ff)summary(ff)@AICtab[3]))
k=sapply(bg.fits,function(mm)length(fixef(mm)))+2
stat.res=data.frame(unlist(bg.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"growth_res_BT.csv")

# could get away with model #33 - veg.bio07 + bolt08 as a growth model 79% of the weight, followed by #25 and #30

# to get the parameters for use in the model
fixef(bg.fits[[33]])
# variances and covariances
vcov(bg.fits[[33]])
# standard errors of the estimates
sqrt(diag(vcov(bg.fits[[33]])))

# check residuals
# plot(fitted(bg.fits[[33]]),resid(bg.fits[[33]]),pch=19+as.numeric(grow.bull[-remove,"bolt08"]))# no obvious non-linearities
# qqnorm(resid(bg.fits[[33]]))
# qqline(resid(bg.fits[[33]]))
# plot(fitted(bg.fits[[33]]),sqrt(resid(bg.fits[[33]])^2),pch=19+as.numeric(grow.bull[-remove,"bolt08"])) # no obvious non-linearities

# make a line plot of predicted population effects
nd = data.frame(veg.bio07=rep(seq(-3.5,1.6,0.1),times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=104),levels=levels(grow.bull$Comp)),
		Herb=factor(rep(rep(c("Ambient","Reduced"),3),each=52)),
		bolt08=factor(0,levels=levels(grow.bull$bolt08)))
nd2 = nd
nd2[,"bolt08"] = 1
nd=rbind(nd,nd2)
rm(nd2)
results = matrix(NA,nrow=nrow(nd),ncol=length(bg.fits))
for (i in 1:length(bg.fits)){
	# generate predictions from a lmer model using the data from inside the model, mostly. 
	fm = bg.fits[[i]]
  results[,i] = predict(fm,newdata=nd,re.form=~0)
}
# model averaged prediction for veg.bio08 as a function of 
nd[,"pveg.bio08"] = apply(results,1,function(x)sum(x*wbic))


# xyplot(log(veg.bio)~veg.bio07|Herb+Comp,data=grow.bull,groups=bolt08,
# 	panel=function(...){
# 		panel.xyplot(...)
# 		pn = which.packet()
# 		pick = nd$Herb == levels(nd$Herb)[pn[1]] & nd$Comp == levels(nd$Comp)[pn[2]] & nd$bolt08 == 0
# 		llines(nd[pick,"veg.bio07"],nd[pick,"pveg.bio08"],col="blue")
# 		pick = nd$Herb == levels(nd$Herb)[pn[1]] & nd$Comp == levels(nd$Comp)[pn[2]] & nd$bolt08 == 1
# 		llines(nd[pick,"veg.bio07"],nd[pick,"pveg.bio08"],col="magenta")
# 
# 	}
# )

coef.names = names(fixef(bg.fits[[1]]))
coef.results = matrix(NA,nrow=length(bg.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(bg.fits),ncol=length(coef.names))
for (i in 1:length(bg.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(bg.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(bg.fits[[i]])))
	names(se) = names(fixef(bg.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}
coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept.nobolt = data.frame(veg.bio07=rep(0,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(grow.bull$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)),
		bolt08=0)
mm.intercept.nobolt = model.matrix(as.formula(bg.models[[1]]),data=nd.intercept.nobolt)

nd.intercept.bolt = data.frame(veg.bio07=rep(0,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(grow.bull$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)),
		bolt08=1)
mm.intercept.bolt = model.matrix(as.formula(bg.models[[1]]),data=nd.intercept.bolt)


nd.slope.nobolt = data.frame(veg.bio07=rep(1,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(grow.bull$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)),
		bolt08=0)
mm.slope.nobolt = model.matrix(as.formula(bg.models[[1]]),data=nd.slope.nobolt)
mm.slope.nobolt[,c(1,3:6,11:15)] = 0
nd.slope.bolt = data.frame(veg.bio07=rep(1,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(grow.bull$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)),
		bolt08=1)
mm.slope.bolt = model.matrix(as.formula(bg.models[[1]]),data=nd.slope.bolt)
mm.slope.bolt[,c(1,3:6,11:15)] = 0


int.results = array(NA,dim=c(length(bg.fits),nrow(nd.intercept.nobolt),2))
int.se = array(NA,dim=c(length(bg.fits),nrow(nd.intercept.nobolt),2))
slope.results = array(NA,dim=c(length(bg.fits),nrow(nd.slope.nobolt),2))
slope.se = array(NA,dim=c(length(bg.fits),nrow(nd.intercept.nobolt),2))

for (i in 1:length(bg.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(bg.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(bg.fits[[i]]))
	int.results[i,,1] = mm.intercept.nobolt %*% fef
	int.results[i,,2] = mm.intercept.bolt %*% fef
	slope.results[i,,1] = mm.slope.nobolt %*% fef
	slope.results[i,,2] = mm.slope.bolt %*% fef
	for (j in 1:nrow(mm.intercept.bolt)){
		pick = which(mm.intercept.nobolt[j,]>0)
		int.se[i,j,1] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope.nobolt[j,]>0)
		slope.se[i,j,1] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.intercept.bolt[j,]>0)
		int.se[i,j,2] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope.bolt[j,]>0)
		slope.se[i,j,2] = sqrt(sum(vcf[pick,pick]))
	}
}

int.avg = apply(int.results,c(2,3),function(x)sum(x*wbic,na.rm=TRUE))
slope.avg = apply(slope.results,c(2,3),function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,c(2,3),int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,c(2,3),function(x)sum(wbic*x,na.rm=TRUE))    ##use this for bootstrap
t1 = sweep(slope.results,c(2,3),slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,c(2,3),function(x)sum(wbic*x,na.rm=TRUE))   ##use this for bootstrap

#residual variance, model averaged
ResVar=sum(sapply(bg.fits,function(fm)attr(VarCorr(fm),"sc")) * wbic)^2

IPM.parameters=data.frame(species=c(rep("BT",6),rep("TT",6)), herb=rep(c(rep("Ambient",3),rep("Reduced",3)),2),comp=rep(c("Low","Medium", "High"),4),
	growth.intercept.nobolt=rep(NA,12), growth.slope.nobolt=rep(NA,12),ResVar.nobolt=rep(NA,12),
growth.intercept.nobolt.se=rep(NA,12),growth.slope.nobolt.se=rep(NA,12),
	growth.intercept.bolt=rep(NA,12), growth.slope.bolt=rep(NA,12),ResVar.bolt=rep(NA,12),
growth.intercept.bolt.int.se=rep(NA,12),growth.slope.bolt.se=rep(NA,12),
  	surv.intercept=rep(NA,12), surv.intercept.se=rep(NA,12),surv.slope=rep(NA,12),surv.slope.se=rep(NA,12),
	bolt.intercept=rep(NA,12), bolt.intercept.se=rep(NA,12),bolt.slope=rep(NA,12),bolt.slope.se=rep(NA,12),
	recruit.ave=rep(NA,12), recruit.se=rep(NA,12),
	size.dist.mean=rep(NA,12), size.dist.unconditional.se=rep(NA,12),size.dist.resse=rep(NA,12),
	zero.seed.intercept=rep(NA,12), zero.seed.intercept.se=rep(NA,12),zero.seed.slope=rep(NA,12),zero.seed.slope.se=rep(NA,12),
 	seed.intercept=rep(NA,12), seed.intercept.se=rep(NA,12),seed.slope=rep(NA,12),seed.slope.se=rep(NA,12)
)

#IPM.parameters$growth.intercept.nobolt[1]=lm(pveg.bio08~ veg.bio07, data = nd[(nd$Comp=="Low") & (nd$Herb=="Ambient")  & (nd$bolt08=="0"),])[[1]][1]
#IPM.parameters$growth.slope.nobolt[1]=lm(pveg.bio08~ veg.bio07, data = nd[(nd$Comp=="Low") & (nd$Herb=="Ambient")  & (nd$bolt08=="0"),])[[1]][2]

for (count in 1:6) {
j=nd.intercept.nobolt$Herb[count]
i=nd.intercept.nobolt$Comp[count]
IPM.parameters$growth.intercept.nobolt[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=int.avg[count,1]
IPM.parameters$growth.intercept.nobolt.se[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=int.avg.se[count,1]

IPM.parameters$growth.intercept.bolt[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=int.avg[count,2]
IPM.parameters$growth.intercept.bolt.se[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=int.avg.se[count,2]

IPM.parameters$growth.slope.nobolt[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=slope.avg[count,1]
IPM.parameters$growth.slope.nobolt.se[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=slope.avg.se[count,1]

IPM.parameters$growth.slope.bolt[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=slope.avg[count,2]
IPM.parameters$growth.slope.bolt.se[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=slope.avg.se[count,2]

IPM.parameters$ResVar.nobolt[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=ResVar
IPM.parameters$ResVar.bolt[IPM.parameters$species=="BT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=ResVar
}


#######################################
# tall thistle biomass as a function of leaves 
#######################################

# make a unique "plot" identifier
grow.tall$plot = paste(grow.tall$block,grow.tall$plot,sep=".")

lgt.mm0 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (1|block) + (1|plot),data=grow.tall)
lgt.mm1 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (1|block),data=grow.tall)
lgt.mm2 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (Herb|block) + (1|plot),data=grow.tall)
# convergence issues when combining comp|plot again
lgt.mm3 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (1|block) + (Comp|plot),data=grow.tall)
lgt.mm4 <- lmer(log(veg.bio)~(log(no.leaves08) + Herb + Comp + bolt08)^2 + (1|plot),data=grow.tall)
sapply(list(lgt.mm0,lgt.mm1,lgt.mm2,lgt.mm3,lgt.mm4),BIC)
# go with block and plot effect - clear winner

# create model set
tb.models = list("~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + 	Herb:Comp + Herb:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + 	Herb:Comp + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + Herb:bolt08 + Comp:bolt08 ",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp + log(no.leaves08):bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Herb + log(no.leaves08):bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + Herb:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp +  log(no.leaves08):Herb + Herb:Comp",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp + log(no.leaves08):bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Comp",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):Herb",
		"~log(no.leaves08) + Herb + Comp + bolt08 + log(no.leaves08):bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + Herb:Comp",
		"~log(no.leaves08) + Herb + Comp + bolt08 + Herb:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + Comp + bolt08",
		"~log(no.leaves08) + Herb + bolt08 +  log(no.leaves08):Herb + log(no.leaves08):bolt08 + Herb:bolt08",
		"~log(no.leaves08) + Herb + Comp + log(no.leaves08):Comp +  log(no.leaves08):Herb + Herb:Comp",
		"~log(no.leaves08) + Comp + bolt08 + log(no.leaves08):Comp + log(no.leaves08):bolt08 + Comp:bolt08",
		"~log(no.leaves08) + Herb + bolt08 +  log(no.leaves08):Herb + log(no.leaves08):bolt08",
		"~log(no.leaves08) + Herb + Comp + log(no.leaves08):Comp +  log(no.leaves08):Herb",
		"~log(no.leaves08) + Comp + bolt08",
		"~log(no.leaves08) + Herb + bolt08",
		"~log(no.leaves08) + Herb + Comp",
		"~log(no.leaves08) + Comp + bolt08 + log(no.leaves08):Comp + log(no.leaves08):bolt08",
		"~log(no.leaves08) + Herb +  log(no.leaves08):Herb",
		"~log(no.leaves08) + Comp + log(no.leaves08):Comp",
		"~log(no.leaves08) + bolt08 + log(no.leaves08):bolt08",
		"~log(no.leaves08) + Herb",
		"~log(no.leaves08) + Comp",
		"~log(no.leaves08) + bolt08",
		"~log(no.leaves08)",
		"~1")

tb.fits = lapply(tb.models,function(ff)lmer(paste("log(veg.bio)",ff,"+ (1|block) + (1|plot)"),data=grow.tall,REML=FALSE))
bic = unlist(sapply(tb.fits,BIC))
dbic = bic-min(bic)
tb.wbic = exp(-dbic/2)/sum(exp(-dbic/2))

data.frame(unlist(tb.models),dbic,tb.wbic)[order(dbic),]
# models 30, 31 tied, then 25, 34, 22, 33, 14, and 28 all do OK.
# # make some descriptive plots 
# library(lattice)
# xyplot(log(veg.bio)~log(no.leaves08)|bolt08+Comp,data=grow.tall,groups=Herb)
# xyplot(log(no.leaves08)~log(no.leaves07)|bolt08+Comp,data=grow.tall,groups=Herb,
# 	panel=function(...){
# 		panel.xyplot(...)
# 		panel.abline(a=0,b=1)
# 	}
# )


# ===================================================================
# growth model using veg.bio.07 - Tall thistle
# ==================================================================
nd = grow.tall
nd$no.leaves08  = nd$no.leaves07
nd$bolt08 = factor(0,levels=levels(grow.tall$bolt08))
results = matrix(NA,nrow=nrow(grow.tall),ncol=length(tb.fits))
for (i in 1:length(tb.fits)){
	# generate predictions from a lmer model using the data from inside the model, mostly. 
	fm = tb.fits[[i]]
	results[,i] = predict(fm,newdata=nd)
}
# model averaged prediction for veg.bio07 parameter
grow.tall[,"veg.bio07"] = apply(results,1,function(x)sum(x*wbic))
# xyplot(log(veg.bio)~veg.bio07|Herb+Comp,data=grow.tall,groups=bolt08,
# 	panel=function(...){
# 		panel.xyplot(...)
# 		panel.abline(a=0,b=1)
# 	}
# )


# create model set
tg.models = list("~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + veg.bio07:bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + veg.bio07:bolt08 + 	Herb:Comp + Herb:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + veg.bio07:bolt08 + 	Herb:Comp + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + veg.bio07:bolt08 + Herb:bolt08 + Comp:bolt08 ",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp + veg.bio07:bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Herb + veg.bio07:bolt08 + 	Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 +  veg.bio07:Herb + veg.bio07:bolt08 + Herb:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp +  veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp + veg.bio07:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + Herb:Comp + Herb:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Comp",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + bolt08 + veg.bio07:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + Herb:Comp",
		"~veg.bio07 + Herb + Comp + bolt08 + Herb:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + Comp + bolt08",
		"~veg.bio07 + Herb + bolt08 +  veg.bio07:Herb + veg.bio07:bolt08 + Herb:bolt08",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp +  veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Comp + bolt08 + veg.bio07:Comp + veg.bio07:bolt08 + Comp:bolt08",
		"~veg.bio07 + Herb + bolt08 +  veg.bio07:Herb + veg.bio07:bolt08",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp +  veg.bio07:Herb",
		"~veg.bio07 + Comp + bolt08 + veg.bio07:Comp + veg.bio07:bolt08",
		"~veg.bio07 + Herb + bolt08",
		"~veg.bio07 + Herb + Comp",
		"~veg.bio07 + Comp + bolt08",
		"~veg.bio07 + Herb +  veg.bio07:Herb",
		"~veg.bio07 + Comp + veg.bio07:Comp",
		"~veg.bio07 + bolt08 + veg.bio07:bolt08",
		"~veg.bio07 + Herb",
		"~veg.bio07 + Comp",
		"~veg.bio07 + bolt08",
		"~veg.bio07",
		"~1")



grt.mm0 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(1|block) + (1|plot),data=grow.tall,REML=TRUE)
grt.mm1 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(1|plot),data=grow.tall,REML=TRUE)
grt.mm2 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(1|block),data=grow.tall,REML=TRUE)
grt.mm3 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(Herb|block) + (1|plot),data=grow.tall,REML=TRUE)
grt.mm4 = lmer(log(veg.bio)~veg.bio07+Herb+Comp+bolt08 + 
		veg.bio07:Herb + veg.bio07:Comp + veg.bio07:bolt08 +
		Herb:Comp +  Comp:bolt08 + Herb:bolt08 + 
		(1|block) + (Comp|plot),data=grow.tall,REML=TRUE)

sapply(list(grt.mm0,grt.mm1,grt.mm2,grt.mm3,grt.mm4),BIC)
# OK, block effect only
tg.fits = lapply(tg.models,function(ff)lmer(paste("log(veg.bio)",ff,"+ (1|block)"),data=grow.tall,REML=FALSE))
bic = sapply(tg.fits,BIC)		
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
loglik = sapply(tg.fits,logLik)
k=sapply(tg.fits,function(mm)length(fixef(mm)))+2
stat.res=data.frame(unlist(tg.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"growth_res_TT.csv")

data.frame(unlist(tg.models),dbic,wbic)[order(dbic),]
# and the winner is model 25 hands down

# to get the parameters for use in the model
fixef(tg.fits[[25]])
# variances and covariances
vcov(tg.fits[[25]])
# standard errors of the estimates
sqrt(diag(vcov(tg.fits[[25]])))

# check residuals
# plot(fitted(tg.fits[[25]]),resid(tg.fits[[25]]),pch=19+as.numeric(grow.tall[,"bolt08"]))# no obvious non-linearities
# qqnorm(resid(tg.fits[[25]]))
# qqline(resid(tg.fits[[25]]))
# plot(fitted(tg.fits[[25]]),sqrt(resid(tg.fits[[25]])^2),pch=19+as.numeric(grow.tall[,"bolt08"])) # no obvious non-linearities
# 
# make a line plot of predicted population effects
nd = data.frame(veg.bio07=rep(seq(-3.5,1.6,0.1),times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=104),levels=levels(grow.tall$Comp)),
		Herb=factor(rep(rep(c("Ambient","Reduced"),3),each=52)),
		bolt08=factor(0,levels=levels(grow.tall$bolt08)))
nd2 = nd
nd2[,"bolt08"] = 1
nd=rbind(nd,nd2)
rm(nd2)
results = matrix(NA,nrow=nrow(nd),ncol=length(tg.fits))
for (i in 1:length(tg.fits)){
	# generate predictions from a lmer model using the data from inside the model, mostly. 
	fm = tg.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~0)
}
# model averaged prediction for veg.bio08 as a function of 
nd[,"pveg.bio08"] = apply(results,1,function(x)sum(x*wbic))

# xyplot(log(veg.bio)~veg.bio07|Herb+Comp,data=grow.tall,groups=bolt08,
# 	panel=function(...){
# 		panel.xyplot(...)
# 		pn = which.packet()
# 		pick = nd$Herb == levels(nd$Herb)[pn[1]] & nd$Comp == levels(nd$Comp)[pn[2]] & nd$bolt08 == 0
# 		llines(nd[pick,"veg.bio07"],nd[pick,"pveg.bio08"],col="blue")
# 		pick = nd$Herb == levels(nd$Herb)[pn[1]] & nd$Comp == levels(nd$Comp)[pn[2]] & nd$bolt08 == 1
# 		llines(nd[pick,"veg.bio07"],nd[pick,"pveg.bio08"],col="magenta")
# 	}
# )


coef.names = names(fixef(tg.fits[[1]]))
coef.results = matrix(NA,nrow=length(tg.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(tg.fits),ncol=length(coef.names))
for (i in 1:length(tg.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(tg.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(tg.fits[[i]])))
	names(se) = names(fixef(tg.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}
coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept.nobolt = data.frame(veg.bio07=rep(0,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(grow.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)),
		bolt08=0)
mm.intercept.nobolt = model.matrix(as.formula(tg.models[[1]]),data=nd.intercept.nobolt)

nd.intercept.bolt = data.frame(veg.bio07=rep(0,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(grow.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)),
		bolt08=1)
mm.intercept.bolt = model.matrix(as.formula(tg.models[[1]]),data=nd.intercept.bolt)


nd.slope.nobolt = data.frame(veg.bio07=rep(1,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(grow.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)),
		bolt08=0)
mm.slope.nobolt = model.matrix(as.formula(tg.models[[1]]),data=nd.slope.nobolt)
mm.slope.nobolt[,c(1,3:6,11:15)] = 0
nd.slope.bolt = data.frame(veg.bio07=rep(1,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(grow.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)),
		bolt08=1)
mm.slope.bolt = model.matrix(as.formula(tg.models[[1]]),data=nd.slope.bolt)
mm.slope.bolt[,c(1,3:6,11:15)] = 0

int.results = array(NA,dim=c(length(tg.fits),nrow(nd.intercept.nobolt),2))
int.se = array(NA,dim=c(length(tg.fits),nrow(nd.intercept.nobolt),2))
slope.results = array(NA,dim=c(length(tg.fits),nrow(nd.slope.nobolt),2))
slope.se = array(NA,dim=c(length(tg.fits),nrow(nd.intercept.nobolt),2))

for (i in 1:length(tg.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(tg.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(tg.fits[[i]]))
	int.results[i,,1] = mm.intercept.nobolt %*% fef
	int.results[i,,2] = mm.intercept.bolt %*% fef
	slope.results[i,,1] = mm.slope.nobolt %*% fef
	slope.results[i,,2] = mm.slope.bolt %*% fef
	for (j in 1:nrow(mm.intercept.bolt)){
		pick = which(mm.intercept.nobolt[j,]>0)
		int.se[i,j,1] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope.nobolt[j,]>0)
		slope.se[i,j,1] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.intercept.bolt[j,]>0)
		int.se[i,j,2] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope.bolt[j,]>0)
		slope.se[i,j,2] = sqrt(sum(vcf[pick,pick]))
	}
}

int.avg = apply(int.results,c(2,3),function(x)sum(x*wbic,na.rm=TRUE))
slope.avg = apply(slope.results,c(2,3),function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,c(2,3),int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,c(2,3),function(x)sum(wbic*x,na.rm=TRUE))
t1 = sweep(slope.results,c(2,3),slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,c(2,3),function(x)sum(wbic*x,na.rm=TRUE))

#residual variance, model averaged
ResVar=sum(sapply(tg.fits,function(fm)attr(VarCorr(fm),"sc")) * wbic)^2



for (count in 1:6) {
j=nd.intercept.nobolt$Herb[count]
i=nd.intercept.nobolt$Comp[count]
IPM.parameters$growth.intercept.nobolt[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=int.avg[count,1]
IPM.parameters$growth.intercept.nobolt.se[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=int.avg.se[count,1]

IPM.parameters$growth.intercept.bolt[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=int.avg[count,2]
IPM.parameters$growth.intercept.bolt.se[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=int.avg.se[count,2]

IPM.parameters$growth.slope.nobolt[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=slope.avg[count,1]
IPM.parameters$growth.slope.nobolt.se[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=slope.avg.se[count,1]

IPM.parameters$growth.slope.bolt[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=slope.avg[count,2]
IPM.parameters$growth.slope.bolt.se[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=slope.avg.se[count,2]

IPM.parameters$ResVar.nobolt[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=ResVar
IPM.parameters$ResVar.bolt[IPM.parameters$species=="TT"&IPM.parameters$comp==i & IPM.parameters$herb==j]=ResVar
}
 
write.csv(IPM.parameters,"IPM.Parameters.csv",row.names = FALSE)